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Non-linear power spectra of dark and luminous 
matter in halo model of structure formation 
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The paper represents a comprehensive treatment of late stages of large-scale structure evolution 
within the framework of halo model. A number of modifications to basic theory are suggested. We 
have engineered simple yet accurate approximation to relate the amplitude of non-linear spherical 
density perturbation to the one of the linear. The theory for final stages of spherical overdensity 
evolution is revised in order to re-evaluate the dependences of collapse and critical overdensity pa- 
rameters, 5 co i, 5 ta and S m i n , on redshift and other cosmological parameters. A new technique is 
proposed for straightforward computation of halo concentration parameter, c, without need to eval- 
uate the z co i. Validity of the technique is proved for a number of ACDM and AWDM cosmologies. 
The parameters for Sheth-Tormen mass function are estimated, as well as new approximation is 
constructed for the dependence of subhalo mass function on initial power spectrum. The modified 
and extended halo model is applied to the determination the non-linear dark matter and galaxy 
power spectra. The semi-analytical estimation of dark matter power spectrum is verified by com- 
parison with data from numerical simulations. Also the predictions for the galaxy power spectra 
are confronted with 'observed' data from PSCz and SDSS galaxy catalogs. Quite good accordance 
is found. 
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I. INTRODUCTION 

A commonly accepted inflationary paradigm states 
that the large-scale structure (LSS) of the Universe is 
formed through evolution of density perturbations driven 
by gravitational instability. Seeded by quantum fluc- 
tuations, the perturbations are initially small in am- 
plitude, they are expanded by inflation afterwards and 
develop during different epochs when medium is domi- 
nated by density of radiation, matter and dark energy 
consequently. A number of manifestations of large-scale 
structure are available for observation: spatial distribu- 
tion of the luminous objects (galaxies, clusters), cosmic 
microwave background (CMB) anisotropy, the features 
in the spectra of high-redshift quasars (Ly Q -clouds), as 
well as the probes via gravitational lensing. The statis- 
tics, e.g. correlations functions or power spectra, can be 
drawn from observational data to test cosmological mod- 
els and scenarios of large-scale structure formation. 

At some point of evolution the growth of perturbations 
of small scales switches to non-linear regime. The linear 
power spectrum is well-understood, for k < 0.1 hMpc -1 
the transfer functions can be readily computed at per- 
cent accuracy for any feasible cosmology. It is not so for 
the small-scale non-linear power spectrum, mostly due 
to the non- linear terms in equations, along with complex 
physics of processes taking place in baryonic component 
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non-linear evolution - power spectra 

(hydrodynamics, radiation transfer, thermal and chemi- 
cal evolution). At final stages the dynamical relaxation 
leads to galaxy formation. The aim of this paper is to 
develop a technique, capable to build a bridge between 
the initial matter power spectrum and observable galaxy 
power spectra inherently non-linear at small scales. The 
treatment is based on the halo model complemented by 
analytical approximations, the results arc tested and ver- 
ified against data of N-body simulations. 

The N-body simulation is a most powerful tool that 
is used to study the structure formation, especially at 
non-linear scales. The numerical simulation technique 
requires input data, usually this is initial distribution of 
the matterparticles, to be generated by some special pro- 
cedure [l|,|2| for given primordial matter power spectrum. 
Also the cosmology should be defined, namely the nature 
of dark components, values of cosmological parameters, 
etc. 

Within standard scenario introduced in [3(, the gravi- 
tational potential of collisionless dark matter shapes the 
wells to which the previously oscillating baryonic matter 
falls in after the decoupling, until the baryonic matter 
power spectrum reaches the amplitude of dark matter 
power spectrum. At some point, a non-linear pertur- 
bation with amplitude exceeding some critical detaches 
from background expansion, reaches a turnaround point, 
and the collapse starts due to self-gravity. The processes 
of violent relaxation take place subsequently, and the 
virial equilibrium has been set in result, the halo of dark 
matter is formed. The baryonic gas starts to gradually 
cool, fragments into clouds and the luminous tracers in- 
side halos are ignited (see 0, 0] for details) . 

Thus, the spatial distribution of luminous matter 



should strongly correlate with halo one. The correlation 
is inferred from large simulations which take into com- 
plete account the baryon physics and particle dynamics 
of dark matter |6|, 0] > as well as from semi-analytic mod- 
els of galaxy formation [8l4l0| . The numerical technique 
requires the considerable computing power, whereas the 
purely analytical is found unreliable. Thus, the "hybrid" 
approach seems to be optimal, combining the analyti- 
cal model of galaxy formation [IJ, [l2[ with dark mat- 
ter "merger trees" extracted from simulations, meanwhile 
the gas dynamics is treated empirically. However, we are 
still far from full understanding of the physics of various 
processes governing the galaxy formation, we can not still 
yield predictions accurate enough. Therefore, the alter- 
native way is to extract the halo and subhalo statistics 
from simulations, to build the appropriate approxima- 
tions and compare the results with observables from large 
galaxy surveys, such as luminosity, color, morphology, 
stellar masses etc. The descriptions used in such tech- 
niques include the conditional luminosity function (CLF) 
|13l - ll5| . conditional mass function (CMF) [Ty| and stellar 
mass to halo mass relation [I(| LL3I • 

The halo model is a cornerstone of modern theory 
of structure formation. It has been proven to be well- 
motivated, comprehensive and provides plausible expla- 
nation for observational data and results of cosmological 
simulations. It is valid for wide range of cosmologies, as 
long as the statistics of primordial density perturbations 
is Gaussian. It encompasses the non-linear stage of evo- 
lution of density perturbations as well as the dynamical 
relaxation processes assuming that all mass is associated 
with gravitationally bound virialized halo. 

It was shown in papers |18l - [20| that given the halo 
statistics, known internal structure of halo and spatial 
distribution of halos, it is possible to evaluate the non- 
linear power spectrum of dark matter perturbations. At 
the next step, the initial power spectrum can be re- 
constructed by applying the halo model to the data of 
N-body simulations. For example, the halo occupation 
function, p(N\M), was used in [2l| to describe the prob- 
ability of finding N galaxies within a halo of mass M, 
to complement the halo model and to enable thus the 
calculation the non-linear galaxy power spectrum. How- 
ever, regardless of overall success of the halo model 1 in 
giving the description of dark matter and galaxy cluster- 
ing, it is still not complete. For instance, only the latest 
enhancements proposed in [23j, |24| ] take into account the 
halo morphology and internal structure. 

The halo model theory is based on the consideration 
of some individual spherical inhomogeneity evolving from 
early initial state to some point of turnaround, moment 
of collapse and formation of virialized halo. The informa- 
tion about such evolution is compressed into few generic 
parameters, namely the amplitude of linear density per- 



turbation at turnaround point, St a , the same amplitude 
at the moment of collapse of central homogeneous core of 
perturbation, 8 co i, and the relative density of virialized 
halo, A vc (in the literature often is denoted byA c ). All 
these parameters depend on the cosmology, primarily on 
the curvature of 3-space and content of dark energy. 

In section 2 we analyze the evolution of a spherical per- 
turbation in different cosmologies, in order to build ana- 
lytical approximations applied henceforth in halo model. 
Also, the new technique is proposed there for computa- 
tion of concentration parameter c for halo with Navarro- 
Frenk- White density profile. 

The section 3 is devoted to the analysis of dark mat- 
ter and galaxy clustering, by means of the halo model 
modified by including statistics drawn from simulations. 
In section 4 the model is applied to estimation of the 
initial power spectrum of dark matter from 'observable' 
galaxy power spectra, which are derived in |25l |26[ from 
PSCz and SDSS galaxy catalogs respectively. The de- 
pendence of subhalo mass function (see 27]) on initial 
power spectrum is investigated, the new approximations 
are proposed for it. The non-linear stage of evolution is 
discussed therein. The conclusions are presented in sec- 
tion 5. All computations were performed for ACDM (cold 
dark matter with A-term) and AWDM (warm dark mat- 
ter with A-term) cosmological models, some bulk math- 
ematical derivations are moved into appendices. 



II. FORMATION OF INDIVIDUAL SPHERICAL 
HALO 

A. Background cosmological model 

We assume that Universe is a homogeneous and 
isotropic space-time with line-element in the general form 

ds 2 = dt 2 -a 2 (t)M 2 {R)dR 2 -a 2 {t) X 2 {R){d6 2 +sm 2 6dip 2 ), 

(1) 
where t is cosmological time and R, 9, ip are spherical co- 
ordinates of 3-space with curvature K (c = 1 here and 
henceforth). The functions x(i?) and M(R) are related 
by 



dx(R) 

dR 



(1-K X 2 (R)) 1/2 M(R). 



Either x(i?) or M(R) can be defined deliberately because 
an arbitrary transformation of the kind R — > R(R) is 
allowed. For example, setting ftA(R) = 1 leads to the 
commonly used line-element with 



R for K = 0; 

X(R) = { Wl{s/KR)/VK for K > 0; 

sinh(y/\K\R)/y/\K\ for K < 



(2) 



1 The details of halo model can be found in 122 



Assuming x(R) = R we obtain other well-known form 
with M{R) = (1 - KR 2 )- 1 ' 2 . 



Another assumption is that expansion of the Universe 
is described by the Einstein's equations with cosmological 
constant, 

Q^v = fcjtu/ - -zpv.v'R- = %-kGT^v + .g M „A, (3) 

where TZ^ is Ricci tensor in the space-time with met- 
ric g^y defined by flj) and T^ v is an energy-momentum 
tensor of matter. After recombination epoch the energy- 
momentum tensor is the one of pressureless ideal fluid 
with energy density p m (t), T>? = diag{p m (t), 0, 0, 0}, 
expressed in comoving reference system. The matter in- 
cludes the baryons and dark matter, which assumed to be 
of the cold or warm type. The dimensionless quantities 
of matter density, curvature and cosmological constant 
are defined as 



h 



« "m — 



3#o 2 



n A 



A 
3Hl> 



-K 

^ = 772": W 



(p^j is mean density of dust matter and H is Hub- 
ble parameter at present moment of time). The Ein- 
stein's equations Q along with the differential energy- 
momentum conservation law lead to the Friedmann equa- 
tions in the commonly used form 



H = H (n m /a 3 + n K /a 2 + ilA 



H (V/« 3 



H 2 \2 



0/ 



(5) 
(6) 



with H = a/a as Hubble parameter and q = —ad/a 2 as 
deceleration parameter at any moment of time (overdot 
denotes d/dt). For current epoch t — t , the equations 
are simple: f2 m + VLk + &A = 1 and qo = ^l m j1 — Qa- 



B. Spherically symmetrical perturbation 

Now let us analyze the non-linear stage of evolution 
of perturbation with spherical symmetry in different cos- 
mologies in order to determine the dependences of the 
characteristic amplitudes on curvature, cosmological con- 
stant and matter density. Within Tolman's approach 
28] the Einstein's equations (J3|) have analytical solu- 
tions for spherically symmetrical inhomogeneity in the 
synchronous gauge (i.e. in the frame comoving to the 
dust-like matter component). The space-time interval in 
this case has the form 



ds i = dt 2 - e^'^dR 2 - r 2 (t, R) (d0 2 + sin 2 



. (7) 



where t is proper time of observers at any R. Henceforth 
it is designated by r after rescaling to units of Hq 1 . 
For the metric ([7]) the Einstein equation Q® = yields 
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FIG. 1: Top panel: The initial profiles of spherical perturba- 
tion of density, S P (R), and mass, Sm(R), both renormalized to 
their central values. The relative gravitational mass defect, 
AMr/Mr, is presented in the units of relative mass defect 
for homogeneous perturbations with the equivalent local cur- 
vature Qf(R): HqR?Q.j(R). Bottom panel: The evolution of 
density profile at linear and weakly non-linear stages, as well 
as for non- linear stage (collapse). The value of initial ampli- 
tude is tuned in a specific way to provide turnaround point 
at Z — 2.5 and complete collapse of central core at z = in 
cosmology with fl m = 0.3 and Q,a = 0.7. 



(here and below (') = d/dR). Because the product 
r /2 e -A j g a f unc ti on of solely comoving radial coordinate 
R, the last equation is rewritten as 



„/2 -A 



= !-/(#), 



(8) 



where function f(R) is relevant to initial conditions for 
perturbation, its density amplitude and profile at some 
early moment of time, as it will be shown below. 

At some early epoch the amplitudes of metric pertur- 
bations are quite small, so the metric ([7]) converges to ((T|) . 
Consequently, in the limit r — ¥ the terms r(r, R) — > 
a(r)x(R) and e A(r ^ )/2 -► a{r)M(R). As the spherical 
perturbation is local overdensity with some characteristic 
scale R c at R ^> R c the terms r(r. R) « cl(t)x(R) and 
e A(r,fi)/2 w a ( T )M(R). 

Let us introduce two variables x(t, R) = r(r, R)/x(R) 
and y(r, R) = e x ^ T ' R ^ 2 / M(R), which are ratios of corre- 
sponding components of metric (J7J and ([TJ . The depen- 
dence of matter density on time and space coordinates in 
the region of perturbation can be presented as 

P (t,R) = p° m x- 2 (T,R)y- 1 (r,R), 



that follows from the differential energy-momentum con- 
servation law, Ttf = 0. At early epoch, the amplitude 
of perturbation is small and the components of metric 
([7]) can be presented as linear expansion around corre- 
sponding values of background metric ((T|). The intro- 



duced variables x(t,R) and y(r,R) as well as relative 
density perturbation 6 p (r, R) , defined as 



S p (t,R) 



Pm(T,R) 



1, 



(9) 



Pm( T ) 

are represented by linear approximations at this epoch, 
x(t,R) ~ a(r) (l - i (r)a(iZ)^ , (10) 

(11) 
(12) 



y{r,R) ~ a (r) II - -a(r)/3(i?) 

<5 p (r,i?) ~ a(-r)A(i?). 



Here, the functions a(-R) and /3(i?) describe the radial 
dependence of spatial components of metric (|7|) and the 
function A(R) = (2a(i?)+ BjM))/3 is the initial profile 
of density perturbation [29|, [30j . 

From Eq. ([8]) and expansions ([T0|) - ([TT|) two relations 
follow: X '{R)P{R) = x'{R)a{R)+x{R)a'{R) &nda{R) = 
3x" 3 (i?) /* A(R) X 2 (R)x'(R)dR, by which the radial de- 
pendences of metrics components can be determined for 
given profile of initial density perturbation [29], [3(| ■ The 
Einstein equation Q® = 8irGp m + A and approximations 
dinjl-dnjl yield the equality 



4D 8 

-vGplA{R)x 2 {R)x'{R) = M (f( R )x(R) - Kx 3 (R))- 

( 13 ) 

With local curvature parameter defined as 

n f (R) = -f(R)/(H 2 X 2 ) (14) 



by integration of Eq. (|13[) we arrive to the expression 
which relate the local curvature to the radial profile of 
initial density perturbation and curvature of background 
space, 



n f ( x ) = -5ft mX - 3 J A( X )x 2 d X + 



{Ik. 



(15) 



This expression along with definition (TUT) allows the 
function f(R), introduced in ©, to be represented via 
profile of initial density perturbation. 

Thus, the spherically symmetrical perturbation is com- 
pletely defined by defining its initial profile A as function 
of x- I n case °f x(R) = R the expression ([T5"j) is like 



R 

n f (R) = -m m R- 3 / A(R)R 2 dR + n K . 

o 



(16) 



The profile of initial density perturbation can be cho- 
sen deliberately, just for simplicity of derivations we put 
it in the form 



(17) 



MO = So 1 - ?£ 



£ = R/R c , where R c is a characteristic scale and 5q is an 
amplitude of perturbation [3l|. Given that, the Eq. CG 
transforms into 



Q f (0 = --n„Ae-« 



n 



K- 



(18) 



For R ^> R c it is readily seen, that the metric (J7J con- 
verges to (QJ, as {lf(R) asymptotically approaches Qk- 
The general equations for variables x(t, R) and y(r, R) 
are derived from Einstein's equations Q\ = A and Q 2 — 
\Q\ = I A as follows 



ftAX 

{Lav ■ 



1 x 2 1 0/ 
2~ + x~' 1 
xy 1 i 2 y 



n f + R- 



a' 



(19) 

a; 2 2 
(20) 



Here and henceforth a derivative with respect to r is 
denoted by overdot. The initial conditions are adia- 
batic and refer to some early moment of time, n = 
|fi ro - 1/2 3/2 (a, < 1): x, = oj(1 - a i a(i?)/3), Vl = 
a,(l - ai0(R)/2i), x { = d^l - 2a l a(i?)/3), y, = a,(l - 
2ai0(R)/3) and %(i?) is given by (JT5J| 

The first integral of (JT5J) is 
A a; 2 = %• Here F(R) = 2GM R (see 



F(i?) 



(i? 3 J tf 2 )x- 1 

:]), the Mjj 



'<iR is the mass of matter contained within 
the sphere of radius R, and the curvature of space is not 
considered. Hence, we arrive to 



X 



fl\X 



n, 



(21) 



The total (gravitational) mass within the sphere of ra- 
dius R is M R = 4tt j { f Pm r 2 e x ' 2 dR = 4np m J*R 2 (l + 
H^R 2 n f {R))-^ 2 dR. The difference AM R = M R - M R 
is called a gravitational mass defect. It can be either neg- 
ative or positive, depending on the sign of the local curva- 
ture, flf. The relative mass perturbation 6 m = 6M R /M R 
in the sphere of radius R, see [32J , is 



6m(t,R) 



r 3 (r,R)J 

a 3 (r) 
x 3 {t,R) 



6 p {T,R)r 2 (T,R)r'(T,R)dR 



1. 



(22) 



Thus, the value pm{t,R) = /9^x _3 (r, R) is a mean 
matter density in the sphere of radius -R, and p(r, R) = 
p^ n x~ 2 (T,R)y~ 1 (T,R) is a matter density at radius R 
from the center of perturbation. The initial profile of 
density perturbation (|17[) and corresponding profile of 
mass perturbation 6m(R) are shown in Fig. QJ,, renor- 
malized to the central values. The relative mass defect 
is presented in the units of relative mass defect of the 
homogeneous perturbations with the same value of local 
curvature £l f (R) at any R: H^R 2 n f (R). 

Within the approach we are developing here, the mass 
M(R) keeps constant value for any R. It means that 



equations (J19I) and (1201 are valid if flf < VLk and Q' f > 
at any R > 0. In this case the first counter-flow of the 
dark matter component arises in the center of perturba- 
tion at moment of collapse 2 , when r(r c ,R < Rh) = 0, the 
Rh is the radius of central homogeneous region (top) . At 
this moment (t c ) the equations ([T9|) -([20 |) become singular 
for < R < Rh and thus unusable. 

The most of profiles of initial spherical overdensities 
relevant to cosmological applications can be described by 
the central amplitude and the characteristic radius (e.g. 
(TTT)) ). In the model with dust-like matter the growth 
of amplitude of central part does not depend on density 
distribution outside of it, and the time dependence of 
profile (and characteristic radius) is determined by the 
central amplitude. The evolution of density perturbation 
profiles is shown in Fig. [TJd, with initial density profile 
given by (|17[) and value of central amplitude taken to 
provide turnaround point at z — 2.5 and finished collapse 
of central part at z = in the cosmological model with 
fl m = 0.3 and il\ — 0.7. 

The density and mass perturbations, as defined by (|12[) 
and (|22[) correspondingly, in the case of small amplitudes 
(^; 1) can be approximated by 



S p (t,R) ~ S(a(T),R) = A(R)D(a(T)) 



(23) 



5 M {r,R) ~ -5{a{r),R) ^ K -ft ({R) D (o(r))(24) 

Here, D(a(r)) is the growing mode of solution for density 
perturbation in the linear theory, 



5 + 2-5 
a 



a a z a 



(25) 



It can be expressed in the integral form as follows 



X-*l 2 {a)d~a, (26) 



D(a) = -Sl m a- l X l l 2 {a) 



where X(a) = A a 2 + t! m ffl^ + Q K [32ll33j. 

The evolution of a spherical inhomogeneity is accom- 
plished in two stages: 1. The expansion stage (i > 0, 
Pm < 0), or linear and weakly non- linear regime; 2. The 
collapse stage (i < 0, qm > 0), or fully non-linear regime. 
They are separated by the moment of turnaround, when 
x = {p M = 0). 

The equations (|2"Tj) and (|20|) can be integrated numer- 
ically for different cosmologies and values of 5q (& n ex- 
ample of such computations is shown in Fig. [1}. It is 
suitable to treat the evolution of non-linear perturbation 
by comparing it with the evolution of some fictional lin- 
ear perturbation just extrapolated into non-linear stage. 
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FIG. 2: The dependence of non-linear amplitude of spherical 
perturbation on the linear one (top panel). The accuracy of 
the approximation (|27|) for some cosmological models (bottom 
panel). 



So, the time dependence of the perturbation amplitude 
is represented in terms of the ratio of amplitude of linear 
perturbation with initial mean overdensity S at some R 
to 5 C oi, the amplitude of perturbation which is to collapse 
at moment t co i . The dependence of non-linear amplitude 
calculated by integration of equations (j2"Tj) and ([2H)l upon 
the ratio S/8 co i is shown in the top panel of Fig. [2] for 
some cosmological models. 

The fit for the dependence is simple, 



lg [1 + s 



M 



lg [1 - 5/5 col ] + A ■ lg 2 [1 - S/S col ] 



B ■ lg 3 [1 - S/S col ] 



(27) 



coefficients A = 0.0903 and B = 0.0074. The errors of 
the approximation lie within 1% range up to moment of 
finished collapse, when 5m — > °° (bottom panel of Fig. [2]). 
The matter overdensity amplitude 
from 5m for any radius R by relation 



5 P , can be estimated 



1 



'Hi 



l + (5-5)jL\n(l + 5 M ) 



(28) 



C. The turnaround point and minimal amplitude 



2 The final stages of collapse are not properly described by 
Eqs. IU9H - H2UH because the dust approximation fails at some high 
density, when the non-radial fluctuations take place and even the 
dark matter collisionless component undergoes the so-called vi- 
olent relaxation before the virial equilibrium is set. 



It was already mentioned, that the moment of 
turnaround is defined by equation x — 0. The value 
of x at this moment, henceforth x ta , is determined from 
the Eq. (|2ip . which at this moment is simplified to 



fi A a; 3 a + flfXta + ^r, 



0. 



(29) 
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FIG. 3: The dependences of the critical amplitudes, S co i, 5ta 
and Smi n , on the matter content Q m , for different values of the 
curvature parameter: a) £Ik — 0, b) Q.k = 0.2, c) Q.k — —0.2. 



The solution for overdensity (0/ < 0) in model with 
A > is 



( 40/ y J i 

Xta = {-xTj cos 3 arccos 



O m ( 3Qa V 

20 A V 0/ y 



The moment of reaching turnaround point, n a , can be 
evaluated by integration of Eq. ([2"Tj) : 



7"to 



d.r 



(^n-T-i+OA^+O/) 1 ^- 



(31) 



The expression in denominator equals zero at X = x ta , 
nevertheless the integral is finite. 

In cosmology with Q\ < or/and O^ < 0, every per- 
turbation which have positive amplitude (Oj<- — 0/ > 0) 
can reach the turnaround point. However, in cosmology 
with Oa > or/and Qk > 0, only perturbation with am- 
plitude exceeding some critical value can do that. The 
condition of existence of turnaround point follows from 
Eq. (|2ip: the perturbation stops the expansion at x ta if 
function x 2 (x) has extremum (minimum) at x ex t in which 

i 2 (x ext ) < 0. The condition is 0/ < -3(O m /2) 2 / 3 A /3 . 
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FIG. 4: Linear amplitudes of perturbation undergoing dif- 
ferent stages at the time z: a) collapse, 8 co i, b) turnaround, 
S ta , and c) asymptotic approach of turnaround, 8 m in- The 
curvature takes the values Q.k = 0, 0.2, —0.2. 



The equality corresponds to the case of x ta = x ex t and 
x = X = (see Eq. (fl"Tj)) ), therefore the point of ex- 
pansion halt is geometrically a point of inflection not 
turnaround. This equation yields the minimal value of 
anrplitude which a linear perturbation ought to acquire in 
order to arrive at the halt but never collapse completely 



, v 3 K 
t(T) "5 — 



3(0 m /2) 2/3 A /3 



D(a(T)). (32) 



Only perturbations with linear amplitude S(t) > 
$min (t) can reach turnaround point followed by collapse 
and formation of virialized object for cosmologically jus- 
tified time. In the Einstein-de Sitter cosmological model, 
0,„ = 1, Oa = Ox = 0, 8 m i n = and the value of 
Sta does not depend on the moment of turnaround and 
equals 



Sta = 



3tt 

T 



2/3 



1.0624. 



However, for the models with non-zero A-constant and/or 
curvature O^, the value of 5t a does depend on the mo- 



ment of turnaround, z ta = z(r ta ). For instance, at 
given z ta > 5, if the influence of A-constant is negligi- 
ble (0 A < fl m (i + z ta ) 3 , ^A < flK(l + z t a) 2 ) the proper 
linear amplitude of perturbation has simple asymptotic 
dependence on z ta '- 



Sta(zta) 



3tt 
T 



1.0624 + 0.2054 




f2 m (l + z ta ) 



In general case the Sta(zta) is evaluated by the sequence 
of equations 01]), §U\), (J^J), (ESJ), §^, ^ and ^. 
The dependences of values 5t a and 8 m i n on the cosmo- 
logical parameters and redshift z are shown in Figs. |3J- 
|4j Note some peculiarities in Fig. [3] at f2 m ~ 0.032 for 
fix = —0.2. These peculiarities indicate the approaching 
the stationary Einstein model in parameter space, where 
the equality 27^^(1 — fl m — Qk) + 40^ = is satis- 
fied. These models are excluded from consideration by 
enforcing finite positive value for Hubble constant H$, 
according to observations. The case of small tt m was 
omitted too, as far as it corresponds to the models with 
bounce, it is beyond this paper scope. 



D. The collapse of overdensity 

The moment of complete collapse of spherical over- 
density, t co ;, is defined by condition x(R) = 0. At this 
moment the matter density p(R) as well as the magni- 
tude of its relative overdensity Sm(R) become infinite, 
the Eqs. (TTT)]) - (l2"Tj) become singular. Meanwhile the mo- 
ment of collapse t co i can be easily estimated, because 
the solution of Eq. (|2ip is symmetrical with respect to 
the turnaround point, r ta - Thus, the moment of collapse 
is evaluated by 



Tcoi = 2rta = 2 



dx 



(n m x- 1 + n A x 2 + n f ) 1 / 2 ' 



(34) 



It is applied usually along with Qf(R;T co i) (Eqs. (IT^l) . 
(IT7)) and (fM]) ) to determine the linear amplitude, S co i — 
So , of perturbation which collapses at given t co i : 

x I \ _ ^ &k - ttf(R]T co i) 

0coi{Tcoi) = ^ D (a(r coi )) . (35) 

The critical value of amplitude of the perturbation, 
which undergoes the collapse at given moment of time 
t co i , is evaluated for any moment of time r according to 



Scoi(t co i\t) = S co i(t co i)D (a(r)) /D (a(r co/ )) 



(36) 



For the cosmology with fix = the formulas (|35D and 
(|36[) are equivalent to the results presented in [34j . In 



the simplest case of Q r , 
amplitude is 

$coi(t co i\t) = - I — 



1 and Hk = the critical 



2/3 



«(r) 

a{r co i) 



1.68647 



1 + z co i 
l+z ' 



at t = t co i it acquires the canonical value. 

In the models with non-zero A-constant and/or the 
curvature parameter fix , the value of S co i depends on the 
collapse moment r co ;. However, for z co \ > 5, when the 
influence of A-constant is negligible (^a <C ^m(l + z co z) 3 , 
f^A <C ^k(1 + z co i) 2 ), the linear amplitude of collapsing 
perturbation is described by the following asymptotic de- 
pendence on z co i: 



Scol(Zcol) 




"|2/3 N 



n 



K 



*'m \ *■ 



1.6865-0.0264 



O m (l + z col )' 



Zcoi) 

(37) 



Generally speaking, for some particular cosmology the 
dependence S co i{z) evaluated at the moment of collapse 
z co i enables to assess the influence of cosmological pa- 
rameters on the dynamics of perturbation development 
during the evolution of Universe. 

In general case the sequence of Eqs. (TT6"1) . (1501) . (|34l) . 
dMl, (ESJ), ®, (EH) and flU} must be applied to ana- 
lyze the influence of cosmological parameters on the value 

Scol(z). 

The dependences of 8 co i on the parameters of cosmo- 
logical model and the redshift of collapse, z, are presented 
in Figs. [3]|3] along with amplitudes St a and S m i n . One can 
see, that S co i is less sensitive to parameters in comparison 
with 5 ta and S min . As a matter of fact, the magnitude of 
Scoi falls within the narrow range 1.55 < 5 co i < 1.75 for 
the wide range of parameters 0.1 < fl m < 1, < fl\ < 1, 
-0.4 < Q K < 0.4. 



E. Virialization and final parameters of individual 
spherical halo 

Note that true singularity in the center of the pertur- 
bation (i.e. halo) is never reached, because the falling 
of particles is never strictly radial, and small-scale inho- 
mogeneities within the cloud induce additional velocities 
of particles, etc. The process of virialization is far from 
trivial, however, after the relaxation has finished and dy- 
namical equilibrium has been established the virial theo- 
rem provides us with simple relations between the kinetic 
energy and gravitational potential. For instance, for a 
spherical relaxed halo the kinetic energy per unit mass is 
determined by 



1 



T mr /m= -{v) v 



1 dU vi , 
—r 

2 dr 



Within ACDM cosmology U vir — U A \ t 



U, 



m\vir 



-H^Q\x^-^—HQrt m /x V i r and the total energy of isolated 
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FIG. 5: The dependence of matter density of virialized spher- 
ical cloud in units of the critical density at the moment of 
collapse, A V c, on Q m in models with fixed Qa (a) and S2k 
(b). 



Note, that the value of x V i r depends on the local curva- 
ture flf and consequently, on the collapse time, t co \\ the 
virial mass density, p V i r = p^x^, depends also on tool- 
It is convenient to represent the virial density in units of 
critical one at a moment of collapse: 



t-^vc. 



Pvir \J~col) 
Pcr\J~col) 



o, 



i-ffo 



.(t co1 )H 2 (t c , 



(41) 



0) 



For the Einstein-de Sitter model (O m = 1, 0^> 
this ratio does not depend on the collapse moment and 
equals A vc = 18tt 2 ~ 178. In the models with A > 
or/and Hk 7^ it takes different values within the 
range 40-180, as shown in Fig. [5] In the wide range of 
parameters 0.05 < Q m < 1, < Q A < 1, -1 < Qk < 
0.95 it can be estimated with percent accuracy by means 
of approximation 



A vc = A-B\og(Q m + C), 



(42) 



where 



A= 110.89 -12.910a 
B= -113.28 



7.92(0 A ) 2 - 0.03/O m , 
6.840a + 7.13(0 A ) 2 , 



C = 0.787 - 0.1150a - 0.093(O A ) 2 . 

The coefficients A. B, C have been determined by the 
fitting to the exact dependence (I41| . 



dark matter cloud is conserved. By comparing the total 
energy at turnaround point, when kinetic energy is zero 
(E to t = Ut a )i with one at virialization epoch (E to t = 
U V i r + T V i r ) we get the equality 



20 A a4> + -—^- = 9 A x 2 ta + — . 



(38) 



Taking into account the Eq. (|29p for turnaround point 
we obtain the cubic equation for x V i r , 



40 A x^ r + 2QfX V ir + 0.,„ = 0, 



(39) 



which for overdensity (flf < 0) in cosmology with A > 
has real root 



2% 

30 A 



1 

cos < - arccos 
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In A = cosmology, as it follows from Eq. 
the strict equality x V i r = x ta /2 is obeyed. For A > 
x v i r < xta/2, but difference is not large. For example, in 
fiducial model with A = 0.7, O m = 0.3, for perturbation 
which collapses at current epoch, the relative difference 
{\xta — x V ir)/ x V ir is ~ 0.1 and diminishes with either A 
decrease or increase of collapse redshift. Therefore, the 
approximation x V i r ~ Xt a /2 is used in the most applica- 
tions. 



III. HALO MERGES, DENSITY PROFILES AND 
CONCENTRATION PARAMETER 



The large-scale structure of Universe seems to be quite 
complicated, and the N-body numerical simulations are 
the only technique that strives to provide the complete 
explanation of it. While the gravitational clustering of 
dark matter is well-understood, many issues are still rel- 
evant to the problem of connection between spatial dis- 
tribution of dark matter and galaxies. Also, only the- 
oretical analysis can bring some useful insights into the 
problem. The halo model has been proven to be most 
fruitful approach for establishing connections between 
the numerical simulation, theoretical developments and 
Ireal observed structure of the Universe. 



(40) 
A. Relaxation processes and halo density profile 

The violent relaxation is the fastest way for halo to 
reach the dynamical equilibrium. This type of relaxation 
is caused by irregular small-scale density inhomogeneities 
residing in protohalo. The scattering of the individual 
particles and gravitationally bound groups of particles 
on such randomly distributed massive inhomogeneities 
brings the velocity distribution to Maxwell-Boltzmann's 
one, independent of mass. The same is true for mass 



density distribution: 
P( r ) = / mn(m,r) = y mno(m) exp < - 

m m ^ 

3$(r)' 



TO$(r) 
kT 



Ps exp ., 



(43) 



Here, kT = ma 2 /3 is a temperature of particles of mass 
to, velocity dispersion a v depends neither on mass nor 
radius of particles, p s = J2 m rnn o( rn ) * s characteristic 
mass density, 



$(r) = G 



M{r) 



dr, M(r) = 4?r / p{r)r l dr. (44) 



Rewriting Eqs. (|4"5)) in the differential form 



dr 



Mp/Ps) 



3GM(r) 



(45) 



and assuming p(r) — p s (r/r s ) -7 , where r s is some char- 
acteristic radius, we get 7 = 2, p s = (of /67rG)r7 2 , 
and eventually the density profile p(r) = p s (r /r s )~ 2 = 
(al/6irG)r~ 2 . This kind of density profile is called 
isothermal, because the temperature of particles with 
fixed mass is constant over the volume of halo. 

The duration of violent relaxation is estimated as 
(Gpvtr) ~ (200Gp cr ) ^3 Gyrs and it does not 

depend on the mass of the gravitational system. This 
explains how the most of collapsed gravitational systems 
managed to already achieve the state of dynamic equi- 
librium in our Universe with age of ~ 13.7 Gyrs. How- 
ever, despite of expectations, the slopes of the density 
profiles of virialized halos are lower at small radius than 
isothermal. The numerical simulations of cluster forma- 
tion, performed by different teams, provide us with the 
best-fit values for 7, typically varying within the range 
1 



< 7 < 1.5 



The most probable reason for this 
is a phenomenon of orbits separation caused by differ- 
ent gravitational action of central irregularities on the 
collisionless particles which move in orbits with different 
eccentricity. The example of such irregularity is a tem- 
porary increase of the density in the center of protohalo, 
taking place at the moment t co i, when collapsing parti- 
cles encounter for the first time. In result, the particles 
with radial orbits lose kinetic energy so the amplitudes of 
their orbits sharply decrease (see [36( for details) . Mean- 
while a kinetic energy of a particle with circular orbit 
increases, thus the eccentricity of the orbit is enhanc- 
ing. It gives rise to a radial dependence and anisotropy 
of particles velocities distribution as well as to velocity 
dispersion [37H39| . 

In this case, the condition of dynamical equilibrium is 
represented by the Jeans equation 



d In p d In a 2 



din r d In r 



2/3 )=- 



GM(r) 



(46) 



where /3 = 1 — of /a 2 is parameter of anisotropy, at and 
ay are respectively the tangential and radial velocity dis- 
persions, and a 2 = a 2 + <r 2 . 

A universal linear relation between the local radial 
density slope, a(r) — dhip/dlnr, and the local veloc- 
ity anisotropy, /3(r), is identified in [40J |. for the set of 
different systems at equilibrium. According to the re- 
lation, as the anisotropy decreases towards the center of 
halo, the slope of its density profile decreases too. Speak- 
ing more generally, the conclusion is that the emergence 
of anisotropy in velocity dispersion is accompanied by a 
change in the slope of density profile. 

Another possible reason for non-isothermal profile is 
so-called slow relaxation. This type of relaxation is 
driven by interactions of close pairs of particles, and re- 
sults in equalizing the temperatures of particles of differ- 
ent masses. Hence, more massive particles have smaller 
velocity dispersions and appear to be more condensed in 
comparison to the less massive particles, in accordance 
with Boltzmann distribution 



n{m, r) = no (to) exp 



m$(r) 
kT 



(47) 



This phenomenon known as a mass separation. The 
manifestations of slow relaxation vary with particle dis- 
tribution over masses, often it is a decrease of 7. In a 
gravitational system containing N particles with equal 
or close masses the duration of slow relaxation is larger 
than duration of violent relaxation by factor of N/ lg(N). 
For a galaxy cluster consisting of N ~ 10 3 similar galax- 
ies or for a galaxy consisting of N ~ 10 10 stars this 
factor is very large, therefore the duration of slow re- 
laxation for such systems exceeds the age of the Uni- 
verse, ~ 13.7 Gyrs. Nevertheless, in a very heteroge- 
neous system containing particles of different masses the 
mass separation can really be a case, especially if violent 
relaxation is neither complete nor sufficient. 

In N-body cosmological simulations the particles are 
usually assigned with the same mass, although due to 
hierarchical clustering of cold dark matter a large num- 
ber of gravitationally bound subsystems emerge which 
can be considered as particles of different masses. The 
distribution of these subsystems, or subhalos, over the 
masses is far from homogeneous and poorly constrained 
in range. If violent relaxation can not provide the effi- 
cient scattering during protohalo collapse, it results in 
smaller velocities of the final subhalos. Thus, the mass 
separation does occur in halos causing the decrease of 
slope of profiles. 

The halo formation is not finished with virialization. 
The process of merging gradually changes the statistics 
of halos, also the mass of a halo increases. The merging 
process is permanent, although the intensity of merging 
decreases during expansion of the Universe. Thus, the to- 
tal mass and radius of a halo are functions of time. Less 
massive halos are trapped from surroundings by central 
massive halo, they appears to have more energy so they 
are located at larger distances from the center of system. 
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Such spatial separation of collisionless particles of differ- 
ent energies is possible only in inhomogeneous gravita- 
tional field which is constant or slowly changing in time. 
This field can be provided by central part of halo, which 
is left almost untouched after virialization. The shaping 
of halos in expanding Universe visually recalls the emer- 
gence of islands when the level of the surrounding water 
is gradually dropping. Due to permanent processes of 
merging, a halo has no sharp boundaries to separate it 
from environment. Thus, the slope of density profile is 
a < 3 at the outskirts of halo. 

As commonly adopted, halo density profile is described 
by simple general expression, p(r) — / o s (r/r s ) _7 (l + 
r/r s ) 7 ~ Q , with constraints on parameter found by [351 ] 
to be 2.5 < a < 3. The characteristic radius r s spec- 
ifies the distance, at which the slope of density profile 
changes. We use universal NFW density profile hence- 
forth, proposed in [4JJ, this is a special case of general 
profile, with values of slopes fixed as 7 = 1 and a = 3: 



p(r) 



(r/r s )(l+r/r s ) 2 



(48) 



For this density profile the total halo mass diverges log- 
arithmically with r. Therefore we must limit the size of 
each halo by some finite radius. 



B. Statistics of halos 

The definition of halo is a matter of convention. Here, 
the mass of halo is defined as a mass of all matter within 
the volume of radius r V i r . The quantity r V i r is defined 
as a radius of sphere the mean internal density of which 
equals 200 times the value of the critical density. The 
mass of the halo is denoted by M2oo- The mass of halo, 
constrained by sphere with mean internal density equal to 
A vc times the critical one is denoted by Ma- Sometimes 
the mass of halo is defined with factor of 180 times the 
critical density, mass is denoted by Miso in such case. 
The index is omitted when the choice of definition is clear 
from context. 

The ratio of radius r V i r , used at definition of halo, to 
the quantity r s is called concentration parameter (or just 
concentration) and denoted as c. Depending on the defi- 
nition of r v i r , the concentration is appropriately indexed 
as C200, CAi ciso- In order to define the parameters of 
halo profile p s and r s , it is quite enough the define the 
mass of halo and provide the value for concentration. 

The concentration for a halo of fixed mass is a stochas- 
tic variable, distributed by log-normal law, 



p(c\m, z)dc ■■ 



1 



27T<7, 



exp 



In [c/c(m, z)] 



2a 2 . 



dine. (49) 



In this case the variance of concentration eventually does 
not depend on the mass of halo (a c — 0.2 — 0.35, see [4^). 
whereas the mean value of concentration is a function of 
mass and redshift. 



The dependence of mean value of concentration upon 
the mass (called just the mass dependence of concentra- 
tion hereafter) is determined either by analysis of nu- 
merical modeling data or inferred by means of analytical 
methods, see [4l|, |43j, |44j . Since the mass dependence fol- 
lows from initial power spectrum of matter, therefore the 
analytical methods seem to be preferable. The changes 
in mass dependence caused by modifications of shape 
or normalization of initial power spectrum can be eas- 
ily taken into account. 

Analytical techniques to study the mass dependences 
of profile parameters are usually based on the treatment 
of [4l| . The simulations indicate on the growth of profile 
specific density, p s , with decrease of halo mass. The cause 
of this was suggested in [4J| , the less massive but highest 
inhomogeneities tend to collapse earlier. It was assumed 
also that specific density of halo, p s , is proportional to 
the matter density of the Universe, given for the moment 
of collapse, i.e. 



Ps — L/S l m Pcr\± + z col) 



(50) 



where the proportionality parameter C is determined 
from modeling. 

The collapse time is defined by "ad-hoc" condition. 
The collapse is assumed to start at the moment of time, 
when a half of virial mass of the halo consists of progen- 
itors each more massive than a fixed fraction / of final 
mass M. With Press-Schechter formalism this condition 
implies 



erfc 



S C oi(z co i\z) - S co i{z\z) 



y/2(<T*(fM\z) 

leading to the equation: 



a 2 (M\z)) 



1 
2' 



(51) 



Scoi(z co i\z) = 5 col (z\z) + C'y/cr 2 (fM\z) - a 2 (M\z) 
~ Scoiizlz) + C'a{fM\z), (52) 

where a 2 {M\z) = a 2 (M)(D(z)/D(0)) 2 , a 2 (M) = 
a 2 (M\z = 0). Here, z is the moment of halo observation, 
C f=a 0.7, the term a 2 (M\z) was neglected in comparison 
with a 2 (fM\z). The parameters C and / are evaluated 
by the data of modeling. From the analysis of simula- 
tion data, it was established (see [4J|), that the value of 
parameter / sw 0.01 does not change substantially, mean- 
while the value of parameter C ex 10 3 is quite sensitive 
to the changes in parameters of background cosmological 
model and/or initial power spectrum. 

Since the parameters C and / have been defined, the 
specific density of halo, p s , can be estimated from the 
Eqs. (f52"|) and (|5"0|) . in order to evaluate other halo pa- 
rameter r s for given the halo mass M. Such treatment 
lacks the explicit analytical expression for dependence 
$coi{z C oi\z), it ought to be computed for each cosmology 
separately. Therefore it had been proposed in [43[ to 
simplify the equation (1521) : 



(53) 



a(fM A \z col ) = 1.686, 
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/ = 0.01. Also they had proposed the simple approxima- 
tion to calculate the concentration ca = K(l + z co i)/(l + 
z), where the value for constant K = 4 is estimated from 
numerical modelling. 

As it was shown in [44]], the abovementioned treat- 
ment is applicable to the estimation of concentration by 
power spectra in CDM cosmology, however it not so for 
Warm Dark Matter (WDM) cosmology. The wrong de- 
pendences follow for the power spectra of warm dark 
matter, which predict the decrease of power at smaller 
scales. The simulations involving initial power spectrum 
of WDM indicate the increase of concentration with in- 
crease of the mass, whereas for the CDM on the contrary 
the concentration decreases with mass growing. 

Since the slope of power spectrum has the same be- 
havior as the slope of dependence of concentration on 
the mass, it was proposed in [44J to replace the Eqs. (|52|) 
and (|53p with the following: 



Total mass of the halo is 



a e ff(M s \z co i)=C (T 1 , 

a eff (M\z)=a(M\z) 



dlna(M) 

d\nM 



(54) 



where C a 
radius r m 



i 28, the mass M s is determined within the 
= 2.17r s , at which the rotational velocity of 



a particle has a maximum for NFW profile. The following 
estimation has been proposed for the concentration: 



ca 



^vc\Zcol) ' *'m\Z) 
L±vc\Z J ' M m yZ co i } 



1/3 



1 + Z C Ql 

1 + z 



(55) 



This approximation has an obvious drawback, the z co \ 
should be known a priori. It is the major obstacle, as 
the collapse time ought to be computed numerically as 
root of non-linear equation. We propose to overcome it 
by altering a few basic assumptions. First, we assume, 
that specific density of halo is determined primarily by 
the collapse of quite homogeneous central region of proto- 
cloud. The rest of the halo is formed afterwards around 
this core by the collapse of outer shells. The boundary 
of the core is defined as the radius, at which the density 
profile changes its slope, so the mass of the core, M c , is 
confined within these edges. The accurate determination 
of this boundary seems to be impossible, therefore we 
assumed the mass of core to correspond to the mass of 
halo with the radius r c — /3r s , where /3 = 0.5. Thus, the 
value of mean internal density within the radius r c cor- 
responds to the density at the moment when dynamical 
equilibrium is established. 

Given this, we can postulate 

4 , 4 o 3 

M c = -irp mr (z co i)r c = -nA vc (z co i)p cr (z co i)(3 r s , (56) 

were p C r(z C oi) is the critical density at the moment z co i. 
On the other hand, integration of the density profile f4"5)l 
yields 



M c = Aitp^^ 



ln(l + ff) 



P 



1 + /3 



(57) 



M = 4irp s r 
= M, 



ln(l 



1 + c 
ln(l + c) - c/(l + c) 



ln(l + /3)- (3/(1 + /?)• 



(58) 



where c is halo concentration. 

At the next step, in order to evaluate the specific den- 
sity of halo, p s , the new condition of collapse should be 
stated, this time for the part of protohalo with mass of 
M c within the sphere of radius r c , at the observation mo- 
ment. Such condition presum ably would bring together 
the speculations from paper [4l| with corrections pro- 
posed in [4J]. The expressions ([54| conforms to the con- 
dition 



Scol(Zcoi\z) -S co i(z\z) 
^2(a 2 (fM c \z)-a 2 (M c \z)) 



= const, 



(59) 



which is finite value, since / < 1 always. To examine the 
dependence of concentration on the slope and amplitude 
of power spectrum as found in [44]], the term a 2 (M) in 
Eq. (|52|) must be taken into account. Moreover, in the 
case of WDM we have one more reason to keep the term 
a 2 (M) in the equation, since the WDM manifests itself 
in a slow change of r.m.s. of a 2 (M) at small values of 
mass. Thus, in general case the previous approximation 
is not valid. Assuming the value of / to be close to unity, 
the Eq. ([52"]) can be approximated as 



5col{Zcol\z) ~ 5 co i(z\z) + C" 
ld 2 a 2 {M c 



da 2 (M c 



dlnMc 



(1-/) 



2 d\nM? 



<i-fY 



1/2 



D(oy 



(60) 



Due to monotonous dependence of a 2 (M) on In M the 
high-order derivatives may be omitted, and 



Scol(z co i\z) ~ 5 co i(z\z) + g 



da 2 (M c ) 



d\nM r 



1/2 



D{z) 

D(oy 



(61) 



where g is a constant, the value of which can be drawn 
from simulations. 

To confront the Eq. (JBTJ) with that of 0, the ap- 
proximation S co i(z col \z) = 5 co i(z co i\z co i)D(z)/D(z co i) ~ 
1.686D(z)/D(z co i) is used, and the term 5 co i(z\z) is 
dropped as z co i ^> z and 8 co i(z co i\z) ^> 5 co i(z\z) for most 
of observed halos. With these assumptions taken into 
account, the Eq. (ffJTj) is rendered to 

^^-a eff (M c ) = (V2g/1.686)-\ 



a eff {M c ) = a(M c ) 



d\na 
dlnM c 



1/2 



(62) 



The comparison of Eq. (|62p with Eq. (|5"4"|) reveals ap- 
parent difference, namely the powers of derivatives, 1/2 
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in (|62|) versus 1 in (1541) . By the way, the values of 
constants v / 2<?/l-686 and C a are not necessary equal, 
whereas the masses M c and M s are defined with spheres 
with different radii, r c — 0.5r s and r max — 2.l7r s re- 
spectively. We propose to favor Eq. (|6Tj) . it is proved to 
be plausible and well-motivated, comes from condition of 
141], whereas the Eq. (|54|) is based mostly on intuition. 

At the next step we have to estimate the param- 
eters of density profile, so that the relation between 
the critical amplitude S co i(z co i\z) and relative density 
of virialized perturbation A vc (z co i) should be found. 
Note, that according to the definition A„ c (z co ;) = 
n m x-? r (H(z co i)/H )- 2 , where (H(z col )/H ) 2 = ft m (l + 



ol y + n K (i + z col y + n A , 



is given by (j3"9")l . 



The £lf(z co i), local curvature within perturbation col- 
lapsing at the moment z co i can be expressed via critical 
amplitude S co i(z co i\z) due to the Eq. 



Qf{Zcol) = ft 



K 



5„ 5 co i(z co i\z) 



-ft 



D(z) 



(63) 



Further, by evaluating the critical amplitude 
5coi{z C oi\z) for given mass M c with formula (j61[). the 
parameter of profile r s as a function of M c can be found, 
on the base of Eqs. (|39[) and (|63|) . above-mentioned 
definitions and Eq. (|35|l . 



/3 I, 1.163ft m 



1/3 



(64) 



The specific density of halo, p s , is evaluated by Eq. (j57|) . 
The estimation of halo mass can be carried out given the 
known density profile and additional condition. However, 
it is more common to evaluate not just r s and p s , but the 
concentration as a function of halo mass. 

Since the ratio of mean density of halo to specific den- 
sity is a function of concentration, 

MI 3 



Phalo 



4nr 3 - 



ln(l 



1 



(65) 



so the concentration is a function of these densities ratio 

181: 



2 Phalo 

3 Ps 



1-1 Phalo 

2.0 p s 



0.387 



(66) 



It is convenient to express the specific and mean densi- 
ties of halo in units of critical density, i.e. A sc = p s /p C r 
and A/ic = phalo/ per, or in units of mean density of mat- 
ter, as A sm = ps/p m and A hm = Phaio/p m correspond- 
ingly. Thus, phalo/ Ps = Ahc/Asc = Ahm/Asm- Accord- 
ing to the condition, used to define the halo radius r V i r , 
one of the values, either A^ c or A^ m , should be constant 
for all halos, meanwhile either A sc or A sm is evaluated 
by formulas: 



(67) 
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n m Hi ,s 3 /3 


<-±sc 

\ 


xl lr m{z)\n(l + P)~p/{l + py 
1 1 /3 3 /3 


-*sm 


x% ir (1 + zf ln(l + 0) - 0/(1 + /?) 



(68) 
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FIG. 6: The dependence of concentration ca on the mass 
M V i r . The squares represent the data for CDM, and the tri- 
angles do for WDM, both are borrowed from modeling of [44j ]. 
Top panel is plotted for current moment of time, 2 = 0, and 
the bottom one is for z — 2. Solid lines represent our calcu- 
lations. 



The halo mass was calculated by Eq. (I58|) . With analyz- 
ing the numerical simulations, we get g « 6.7. 

The results of calculations of concentration ca are rep- 
resented in Fig.lHl as dependence on the mass M V i r . The 
data of simulation by [44| are plotted for comparison, 
for cosmology with il m = 0.3, Ha = 0.7, h = 0.65 and 
cr 8 = 0.9, both for CDM and WDM. The dependence 
is presented for z — and z — 2. As it is shown, the 
concentration grows with the increase of z co [. Thus, for 
older halos concentration is larger than for younger ones. 
For the ACDM cosmology the halos of smaller masses 
are formed in first turn, subsequently they have larger 
concentration. For the AWDM cosmology the halos with 
smaller masses can be formed later, thus they have the 
smaller value of concentration. 

The evaluated dependence of halo concentration, C200J 
on the mass, M200, is represented by solid lines in 
Fig- El The triangles represent the data from paper [45j . 
taken for cosmological model with parameters provided 
in [HI by 5-year data release of WMAP: Q A = 0.721, 
ft m = 0.279, Q b = 0.0441, h = 0.719, a 8 = 0.796 and 
n s = 0.963. Three plots represent the data calculated for 
the values of redshifts varying z = 0, 1, 2. As it can be 
seen, our calculations agree quite well with the data of 
modeling at all redshifts. 



IV. MASS FUNCTIONS AND POWER 
SPECTRA OF HALOS 

An analytical approach for mass function of clusters of 
galaxies had been suggested in pioneering paper of [47j, 
which is based on the model of spherical collapse and 
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FIG. 7: The dependence of concentration, C200, on the mass 
M2oo- The triangles represent the modelling by [45||. The 
solid lines represent our results. The plots are given for the 
redshifts z = 0, 1, 2 in downward order. 



S to be e„y, = (a(m)/5)/V5, a simple relation was ob- 
tained in [52j to connect ellipsoidal and spherical collapse 
thresholds 



S ec (m,t col ) = 5 col (t co i) 1 + 0.47 



a(m) 



Scol(tcol) 



r-> (70) 

Also, the excursion set model was used in [48[ to es- 
timate the mass function associated with ellipsoidal col- 
lapse 



uF(y) = A{p) (1 + v -*) J — exp {-v/2} , (71) 

where as before v = (8 co i(t co i) / a(m)) 2 , a parameter p ~ 
0.3 and function A(p) = [1 + 2-*T(l/2 - p)/^" 1 ^ 
0.3222 are determined by requirement that all mass is 
gathered within halos, i.e. the integration of F(y) over v 
yields unity. In order to match the data of GIF numerical 
simulations the mass function (|7ip has been parameter- 
ized in [5l| as 



identifies the halos with the peaks of an initially Gaus- 
sian field of density perturbations. This so-called Press- 
Schechter formalism produces the halo mass function to 
describe the distribution of halos over their masses. The 
approach was refined and extended afterwards in [48l450l | 
to allow for the merger histories of dark matter halos. 
The process of halo merging is assumed to be hierarchi- 
cal at the large scales, and described with characteristic 
collapsing mass scale m(t co {) with r.m.s. of density per- 
turbations a(m) — 5 co i(t co i). This mass grows with time 
through merging of halos and in distant future it will 
approach asymptotically the limit value m^ at which 
ff(raoo) = S m i n . For cosmology with 5 m in = the clus- 
tering of dark matter has no end in meaning that all halos 
of the Universe will merge in far future. 



A. Halo mass function 

Authors of [48| have argued that Press-Schechter mass 
function n(m,z), the number density of gravitationally 
bound objects with masses m at redshift z, satisfies the 
condition 



vF{v) = A(p) (1 + (qu)- 



vF(y) 



m 2 n(m, z) dlnm 
~p m din z/ 



— exp{-i//2}, (69) 

Z7T 



where v = (S co i(t co i)/a(m)) 2 and ~p m is the background 
matter density. 

The mass function of Press-Schechter is qualitatively 
correct however it disagrees with the data of N-body sim- 
ulations in some details. Treating the collapsing pertur- 
bations as ellipsoidal rather than spherical diminishes the 
discrepancies (see [511] V Indeed, by assuming the average 
ellipticity of perturbation with mass m and amplitude 



|_exp{-^/2}. (72) 



The additional parameter q was evaluated as q — 0.707, 
later it was re-determined in [53| to be q = 0.75. In 
the framework of excursion set approach the ellipsoidal 
threshold associated with such mass function is as fol- 
lows, see [52J 



6 eq (m,t co i) = q^S co i(tcoi) 1 + 0.5 



a(m) 



q 2 8coi{tcoi) 



1.2N 



(73) 

Two different algorithms are commonly used to iden- 
tify the dark matter halos within data of numerical N- 
body simulations: the friend-of-friend (FOF) algorithm 
[54j and the spherical overdensity (SO) finder 55]. The 
FOF procedure depends on just one free parameter, b, 
which defines the linking length as bn~ 1 ' 3 , where n is 
the average particle density. Thus, in the limit of very 
large number of particles per halo, FOF approximately 
selects the halo as matter enclosed by an isodensity sur- 
face at which p = ~p/b 3 . SO algorithm finds the values of 
the average halo density in spherical volumes of various 
sizes. The criterion for identifying the halo is equality of 
the average density in the sphere to certain value Kp m , 
where ~p m is mean density of matter in sample and k is 
a parameter of algorithm. For the NFW density profile 
these algorithms appear not to be identical as they lead 
to different dependences of the halo concentration c on 
halo mass. Nevertheless, very similar halo mass functions 
have been found in |56|. using SO (k = 180) and FOF 
(b = 0.2) halo finders. 

We refer halo as a gravitationally bound system which 
has reached the dynamical equilibrium, meanwhile both 
SO and FOF finders select the groups of close particles 
regardless of their dynamical properties, halos are treated 
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FIG. 8: Halo mass function from N-body numerical sim- 
ulation by GIF/Virgo collaboration [59( for different red- 
shifts z = 0, 1.05, 1.94, 2.97, 4.04. The dashed lines show the 
Sheth-Tormen approximation with parameters p = 0.3 and 
q = 0.707 [5j|, the dotted lines show the same approximation 
with modified parameters p = 0.32 and q — 0.76. 



simply as dense regions. To divide such halos into viri- 
alized (relaxed) and non-virialized (not relaxed) parts, it 
was suggested in [57} to assess the dynamical state of each 
halo processed by FOF algorithm by means of three ob- 
jective criteria: 1) The substructure mass function f su b', 
2) The center of mass displacement s = \r c — r cm \/r V i r , 
and 3) The virial ratio 2T/U. In [Hj| the r.m.s. of the 
NFW fit to the density profile, p rms , have been utilized 
too. 

As far as virial density p V i r = p cr A vc = ~p m A vm , it 
seems appropriate to use SO halo-finder with k = A vm = 
A vc /£l m . The equality K ~ 180 is valid for any redshift 
in flat ft m — 1 cosmology (this is close to the b ~ 0.2 for 
FOF algorithm), meanwhile for ACDM cosmology with 
Q m = 0.3, 57a = 0.7 the quantities k and b depend on 
redshift: K ~ 97/0.3 ~ 324 (6 ~ 0.164) at z = and 
slowly decrease (increase) to the limit k ~ 180 (b ~ 0.2) 
at high z. However, as it is shown in [56| the shape of 
mass function is invariant if we simply identify clusters 
with a constant linking length, b = 0.2, for all redshift 
and cosmologies. 

The halo mass function derived from N-body simula- 
tions by GIF/Virgo collaboration is plotted in Fig. [5] 
The halo catalogs are built from simulation and made 
available 3 . For each halo, as detected by FOF-algorithm 
(b = 0.2), the catalogs include the mass M200, con- 
fined within the central part of halo with overdensity 



A vm = 200 (see [59| for details). The mass rescaling 
affects slightly the 'observed' mass function. We have 
re-determined the parameters of Sheth-Tormen approxi- 
mation to be p — 0.32 and q = 0.76. As it follows from 
Fig-Hl the refined parameters provide a better fit for data 
than ones from |51j . p — 0.3 and q = 0.707. The mass 
function is defined here as a number density of halos with 
masses exceeding the specified mass m, 



N(> m) — I n{m! ' , z)dm' = 



3 http://www.mpa-garching.mpg.de/GIF 



p m _, .dlnv , 

^-vF(v)——-dm'. 
m dm 

(74) 
Note, that changes in FOF or SO halo finder parameters 
modify also the total number of detected halos. However, 
the mass rescaling influences the mass function shape but 
not the total number of halos. 

According to Press-Schechter formalism [47] halos are 
formed within regions with initial overdensities 5 > S co i 
because these have already reached a dynamical equi- 
librium through the collapse. However, for some re- 
gion to reach overdensity A„ m , the initial overdensity 
can be lower than S co i. Given the non-linear overden- 
sity as 8m — 180 in (|27| . the corresponding initial over- 
density in the units of critical density is calculated as 
q? = S/S co i ~ 0.95. Therefore, the applied SO-algorithm 
(k = A vm ) marks also as halos non-virialized dense re- 
gions formed by growth of perturbations with lower ini- 
tial overdensities, determined by q^5 co i under assump- 
tion of spherical symmetry, the value q < 1 and depends 
on A vm according to (|27|) . It seems reasonable to assume 
that the elliptical 'q-threshold' S eq is directly connected 
to spherical 'q-threshold' qz8 col through ([75]) just as the 
elliptical collapse threshold 8 ec is related to the spher- 
ical collapse threshold S co i through ([70)1 . However the 
value q ~ 0.95 2 ~ 0.90 obtained here for A vm = 180 
differs from the value q ~ 0.75 derived from numerical 
simulations using FOF (b = 0.2) and SO (k = 180) al- 
gorithms. It is known that large halos are mostly spher- 
ical, and hence their mass distribution should be com- 
patible with Press-Schechter mass function. Indeed, if 
the spherical 'q' barrier 0.95<5 CO / (0.9^ for v) is used in 
([69]) instead of 8 co i, a good match is attained to the 
results of numerical simulations for large v and hence 
large masses. However, the Press-Schechter mass func- 
tion predicts the larger number of halos with smaller than 
expected masses. The reason is that less massive halos 
are more elliptical, and therefore according to ([70[) need 
larger initial amplitude to form. 



B. Subhalo mass function 

As it was noted, the clustering of dark matter is a kind 
of hierarchical process, i.e. the halos of small mass have 
come out first, part of them subsequently are merged to 
form massive halos. After initial merging the halo pro- 
genitors gradually scattered during slow relaxation, and 
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the most massive even survived within virialized halo as 
subhalos. A large number of smaller progenitors usu- 
ally mold a smooth part of halo, i.e. diffuse dark mat- 
ter which is not involved into any substructure. The 
considerations on the dark matter halo substructure are 
important to understand the cluste ring of dark matter 
on a small scales ~ 10 3 -=- 10 6 pc J24J, internal struc- 
ture of galaxies [60] and galaxy clusters (see |15Tll7l l61|'). 
The subhalos may also reveal substructures, consequently 
sub-subhalos and so on, the algorithm to identify sub- 
structures through all levels of this hierarchy was pro- 
posed in 27], and applied to determine the mass function 
from numerical simulations. 

The conditional mass function, proposed in [48j, pro- 
vides the description for the merger history of dark mat- 
ter halos, it can be simply derived from (|72[) by putting 
v = [(Scoi(Scoi) - S co i(z co i))/(a(m) - cr(Af))] 2 , where M 
is observed halo mass at z co i and m < M is a mass of 
its progenitor found at some moment z co i > z co i- To 
adapt this finder for the substructure we assume that 
Scoi(zcoi) — o~(m rs ), where m rs is a sort of resolution 
scale, a limitation up to which the substructure is re- 
solved. So, substructures within subhalos of mass below 
m rs are not resolved. 

To provide best-fit to numerical estimation of subhalo 
mass function [271 ] . we multiply our conditional mass 
function by factor of 0.3. It means that halo mass consists 
only by 30% of subhalos and rest 70% is a diffuse dark 
matter which not to be associated with any substructure. 
Hence, we arrive at 



dN . ,d In a 2 (m) 

-= Q.3MvF(u)- 



dln 



dm 



(75) 



wher 



(1.5(5 co/ )7(cr 2 (TO) - cr 2 (Af)), and din a 2 (m) 



is used instead of term din [c 2 (m) — <t 2 (M)] because it 
provides a better accordance when m is close to halo mass 
M. Note that all estimations of mass from numerical 
data were carried out for virial overdensity regions with 

The subhalo mass function depends on the power spec- 
trum via a(m) and a(M), here the type of dark matter is 
involved. For warm and hot dark matter the number of 
low-mass subhalos is expected to be lower in comparison 
to cold dark matter [62J. The previous approximations 
proposed for subhalo mass function in [23, [63[ do not 
include such important dependence. 



C. The dark matter power spectrum 

A process of formation of galaxy as a luminous object 
consists mainly in baryon matter which undergoes con- 
densing and cooling within potential wells formed by in- 
homogeneities of dark matter. Therefore, the observable 
spatial distribution of galaxies should strongly correlate 
with distribution of dark matter. Since the content of 
dark matter (s=s25%) substantially overwhelms the con- 
tent of baryon matter («5%), so the formation of large- 




FIG. 9: Subhalo 
ferent halo masse; 
log 10 (M/M Q ) = 
12.25 (squares). 
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mass function from 
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[2J for three dif- 
14.25 (triangles), 



13.75 (diamonds) and log 10 (M/M e ) = 
Subhalo mass functions for halo masses 
log lo (M/M ) = 14.25 and log 10 (M/M Q ) - 12.25 are rescaled 
by 10 up and down respectively. Solid line show our approx- 
imation (|75[) whereas dashed line show the original approxi- 
mation from 12711 . 



scale structure of Universe is defined primarily by the 
development of dark matter inhomogeneitics. 

Given that, we need to calculate the characteristics of 
distribution of dark matter in order to reconstruct the 
distribution of all matter, i.e. the characteristics of ob- 
servable distribution of galaxies and temperature fluctu- 
ations of cosmic microwave background. 



1. Two-point correlation function and power spectrum, of 
discrete and continuous distributions of matter 

The two-point correlation function provides compre- 
hensive statistical description for inhomogeneous spatial 
distribution of matter. The power spectrum is actually a 
Fourier transformation of the two-point correlation func- 
tion, and can be directly drawn from the spatial distribu- 
tion of matter by means of Fourier-expansion of relative 
fluctuations of density. In the case of continuous distri- 
bution of matter the expansion is following: 

S (r) = P^_~P = (2 7r ) 3 /V 1 / 2 [Sj:e-^d 3 k 



(2^2 



ikf 



(76) 



here, V denotes some volume of periodicity, the result 
does not depend on this quantity. The coefficient of ex- 
pansion can be expressed from (|76|) in the following man- 
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(77) 



The squared Fourier amplitude averaged over the dif- 
ferent directions of vector k is actually what is called the 
power spectrum, V(k) = (|i5j| 2 ). The two-point correla- 
tion function can be derived from known power spectrum, 



t(r) = {6(?)8(f + r f )) = 



(2nf 
V 



£<l%l 2 K fc> 



d Z k(\S % \ 2 )^ = Air \ k 2 dkV(kf^^- ,(78) 



as well as variance of the amplitude within the spherical 
volume of R radius, 

a 2 (R) = <S 2 R >=4n f k 2 V(k)W 2 (kR)dk 

= f A 2 (k)W 2 (kR)d\nk, (79) 

where W(x) — 3(sin(x) — xcos(x))/a; 3 is a window func- 
tion for spherical volume. 

The quantity A 2 (fc) = 47rfc 3, P(fc) is a dimensionless 
power spectrum. The term "dimensionless" is often omit- 
ted, as it is clearly understood from the context of us- 
age. For given correlation function, the power spectrum 
is evaluated by the following formula 



V(k) = 



1 



(2tt) s 



d 3 re 



-iki 



*(r). 



(80) 



The galaxy catalogs (or the data of numerical simula- 
tions) involve discrete distributions of galaxies (or par- 
ticles). In this case the equations (f76| and (f77f should 
be recasted with p(r) = ^miSoif — fj), where TOj is 

i 

the mass of «-th galaxy (particle), 8n(r — ?i) is three- 
dimensional Dirac function. The Eq. (l77|) takes the form 
of 



Sr. 



1 



{2irf/ 2 {m)nV 1 / 



\ii 2_^ 



irijC. 
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here n is spatially averaged number density of galaxies 
or particles, (to) = ^mi/(nV) is the mean mass. 

The power spectrum is related to the correlation func- 
tion by relation [64j In such case 



V(k) 



(to 2 ) 



1 



(27r) 3 n (m) 



(2tt) 



d A re- lKr i[r), (82) 



with (to 2 ) = ^2 m 2 /(nV). The first term in right-hand 
side of (|82")l is a shot noise, denoted by V s hot- It cor- 
responds to the case of discrete distribution, when the 
number density of galaxies (particles) n takes large but 



finite value. At n — > oo, i.e. for the continuous distribu- 
tion the Eqs. ([82]) and (f80| converge. The second term 
in right-hand side of (f82"j) is denoted henceforth as V^{k) 
to emphasize the non-random (correlated) nature of dis- 
tribution. Thus, the Eq. (|82[) can be written in more 
compact form, as V(k) — Vshot + V^{k). 



2. Non-lmear power spectrum in halo model 

Within halo model the distribution of matter is treated 
in a mixed, discrete-continuous manner. The matter 
is segregated into spatially separated halos of different 
mass, somehow distributed over the space, and the mat- 
ter distribution within each halo is described by contin- 
uous density profile (|48j) . Therefore, the power spectrum 
is split into two terms, one to describe the distribution of 
halos, and the second to describe the distribution of mat- 
ter within individual halo. The splitting can be derived 
rigorously taking into account that Fourier-amplitudes of 
density perturbations are (see for details Appendix [SJ 



Sr, = — I m ■ n{m)Sr(m)y(m, k)dm. 
o 



(83) 



Here ~p = p„/a 3 is an average matter density at any 
moment of time determined by scale factor a = (1 + z) -1 
and n(m) is the number density of halo with mass to 
drawn from (|T2[) in comoving coordinates. 

The function y(m, k) is a Fourier-transformation of 
density profile (|48|) expressed by explicit analytical form 



y(m,k) = 



47rp s r 3 
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(84) 



where c is the halo concentration, p s and r s are parame- 
ters of density profile, Si(x) and Ci(x) are integral sine 
and cosine respectively. The halo profile depends on the 
physical coordinates, while the power spectrum is asso- 
ciated with comoving coordinates R — r/a. The term 
47rp s r 3 /m can be expressed via halo concentration pa- 
rameter c using the Eq. (1581) . 

The power spectrum V(k\m, to') of spatial distribution 
of halos with given masses to and to' , is following 



V(k\m, m') 



2^* 
S 



- (^(m)^(m') + S u (m)Si{m') 
V^{k\m, m'). 



(2ir) 3 n(m) 



17 



where 5 m .m' is Kronecker symbol, V^(k\m, m') is Fourier- 
image of two-point cross-correlation function of halos 
with masses m and ml (the angle brackets in the left-hand 
side denote the averaging over the directions of vector k). 
After the series of mathematical transformations, we 



arrive to 



V{k) = 



(2^) 3 7j 



3t;2 



m • n(m)\y(m, k)\ dm 



TABLE I: The parameters of Large Box (LB) and GIF2 sim- 
ulations, available at Max Planck Institute for Astrophysics 
in Garching (http://www.mpa-garching.mpg.de). 



Simul. 


Npar 


L (Mpc/h) 


nip (Msun/h) 


hoft (Kpc/h) 


GIF2 
LB 


400 3 
512 3 


110.0 
479.0 


1.73xl0 9 
6.86xl0 10 


6.6 

30 



_, . m ■ n(m)y(m,k)dm 

P 

o 



(85) 



x / m! ■ n{m!)y{m! ,k)dm''P^(k\m,m'). 


The quantity n(m) is a number density of halos having 
masses m. 

With assumption of linearity of perturbations cross- 
correlation power spectrum can be represented as 
"p£ (k\m, rnf) m b(m)b(m') r Pn n (k), the Vu n (k) is the linear 
power spectrum of spatial distribution of matter, 6(771) is 
the biasing parameter which characterize the skew be- 
tween distributions of halos of mass m and the matter 
distribution. 

Since the distribution of matter is homogeneous at very 
large scales, the expression (|85l) is expected to asymptot- 
ically approach zero for small wave numbers k. However, 
the first term in right-hand side of (|85[) never diminishes. 
The reason is that binning of matter into separate ha- 
los introduces the noise into the procedure, as inherent 
for discrete distribution. The expression for noise is de- 
rived from the first term at (|85|) by letting the distribu- 
tion of halo matter to be homogeneous and substituting 
of Fourier-image of profile y(k, m) by window function 
W{kR), where R = (3m/(47r77 m )) 1/3 . After noise elim- 
ination, the final expression for the power spectrum of 
spatial distribution of dark matter is following: 



V(k) 



1 



— ) n(m) [\y{m, k)\ 2 - W 2 (kR)] dm 



(27r) 3 7 \p 



—61 (rn)n(rn)y(m, k)dm 
P 



Vi m (k). (86) 



We use this equation to calculate the power spectrum of 
dark matter in the wide range of scales and compare it 
with data from N-body numerical simulations. 

Also, the calculations were carried out to confront our 
estimations with results of Large Box and GIF2 N-body 
simulations available at Max Planck Institute for Astro- 
physics in Garching. The simulation team provides data 
files with coordinates, velocities and identification num- 
bers of particles. All particles have equal mass. The 
parameters of simulations (total number of particles, box 
size, mass of particles and smoothing scale) are presented 
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FIG. 10: Dark matter power spectrum obtained from LB (as- 
terisk) and GIF2 (diamonds) simulations. Solid lines repre- 
sent the primordial linear power spectrum from |6g| (lower 
line) and our predictions for non-linear halo model power 
spectrum (upper line). Dash-dotted lines show the 'halo-halo' 
and the 'shot noise' components of the halo model power 
spectrum. Dotted and dashed lines show PD96 [6g] and 
HALOFIT [13 approximations. The ACDM cosmology is 
used with (fi m , fi A , h, a 8 , n s )=(0.3, 0.7, 0.7, 0.9, 1). 



in the Tab. [D The smoothing scale is used to avoid the 
numerical singularities which arise in consequence of par- 
ticles proximity, when the numerical errors is difficult to 
control. 

The higher number density of particles is necessary in 
order to reproduce the structure at small scales. On the 
other hand, the larger volume is necessary in order to 
reproduce the structure on large scales. A pursuit to sat- 
isfy both requirements simultaneously leads to the huge 
numbers of particles, and consequently to the enormous 
amount of computational efforts. So, usually the sepa- 
rate modelings are carried out, to reproduce the matter 
distributions at different scales. In our case the data 
from Large Box simulations cover large volumes, provid- 
ing thus good reproduction at larges scales, whereas the 
data from GIF2 simulations is applicable at small scales, 
having larger average number of particles. 

To evaluate the magnitude of power spectrum at given 
fixed value of wave number k for different fixed direc- 
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FIG. 11: The difference (percentage) between non-linear 
power spectra of ACDM and AWDM models with (Q m , £Ia, 
fib, h, n a , cr 8 )=(0.2711, 0.7289, 0.0451, 0.703, 0.966, 0.809) 
and two values of WDM particle mass: m m dm = 1 keV and 
m w dm = 0.5 keV. 



tions of vector k, denoted by index I, we have estimated 

the sums of Sri = Si cos (^) an d Sii = X)i sm (^)) 
and calculated the power spectra as Vi(k) = (Srf + 
Sif)/((2ir) 3 n 2 V). The net vector was estimated as av- 
erage over all given directions: V{k) = j?Yli=i at^M&)> 
where N is total number of given directions. To eliminate 
the noise, we have precalculated the power spectrum of 
homogeneous distribution of particles, and substracted 
this contribution from the power spectrum for the gen- 
eral case of inhomogeneous distribution of particles. 

The results of comparison are shown in Fig. [TO] As 
it can be seen, the dark matter power spectrum pre- 
dicted by halo model, Eq. flSSJ}, matches the LB/GIF2 
non-linear power spectrum at all scales up to k ~ 100 
h _1 Mpc. Also, the power spectra are shown, calculated 
by PD96 d|| and HALOFIT [67| approximations based 
on the previous halo model, of Hamilton et al. [68| , scal- 
ing relations and fits to numerical simulations. They pro- 
vide a good fit the LB/GIF2 non-linear power spectrum 
at all scales as well. In all cases the analytical approxi- 
mation for linear power spectrum from [65( (lower solid 
line in Fig.lT0"jl is used, normalized to as = 0.9. 

Apparently, the non-linear corrections are not essen- 
tial at k < 0.2, the power spectrum appears to be linear 
therein. At smaller scales the non-linear clustering gains 
domination (dash-dotted line in Fig. [TO]) , and exceeds 
by amplitude the linear power spectrum. Both our ap- 
proximation (|86|l and HALOFIT reproduce such behav- 
ior equally good. Consistency of our power spectrum esti- 
mation to numerical simulations data and HALOFIT ap- 
proximation proves the plausibility of our approach, and, 
on the other hand, gives a deep insight into the process 
of dark matter clustering and foundations of HALOFIT 
approximation, which is widely used. 

To verify the applicability of our halo model in case 
of the WDM, we have confronted the non-linear power 
spectra estimations with the results of numerical simu- 
lations from 16 91 . The initial distribution of warm dark 



matter particles was generated by following linear power 
spectrum of density perturbations 
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(87) 



where 4 v = 1.12. The parameter a depends on the mass 
of WDM particles, m w d m , their density, £l w dm and Hub- 
ble parameter as 



a{m w dm) = 0.049 



IkeV 



1.11 
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>dm 



0.25 



0.11 
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07 



1.22 



(see also [70| and |71|). 

We have computed the non-linear power spectrum of 
matter density perturbations ([56]) at Z = 0.5 for ACDM 
and AWDM cosmologies with (f2 m , ft a, fib, h, n s , 
cr 8 )=(0.2711, 0.7289, 0.0451, 0.703, 0.966, 0.809) and two 
different masses of warm dark matter particles, m w d m , as 
1 and 0.5 keV, the same values for which the N-body sim- 
ulations have been accomplished in [69|]. The Fig.lTTlrep- 
resents the relative discrepancies (in percents) between 
non-linear power spectra of cold and warm dark matter 
(solid lines), the corresponding spectra calculated from 
simulations are plotted along (Fig. 7 in |60|)- Our results 
reveal qualitative consistency with simulations however 
quantitative differences are still noticeable 5 . 

Note the discrepancy in case of warm dark matter halo 
model 72] . depicted by green line in the bottom panels of 
Figure 7 in paper [69J , are also suppressed in comparison 
with numerical simulation and seems to be closer to our 
results. It implies that for WDM models new indepen- 
dent numerical simulations are required to obtain more 
robust statistics as well as the case should be re-analyzed 
in the framework of halo model. Indeed, the amplitude 
of power spectrum is very sensitive to dependence of halo 
concentration parameter on mass of halos, Eqs. (|84|) and 
(|86[) . unfortunately this kind of sensitivity can be veri- 
fied only by high- resolution numerical models (see also 
discussions in I4J, U2 



D. The galaxy power spectrum 

As the baryon gas falls into potential wells of virial- 
ized dark matter halos and subhalos it has been heated 
to virial temperature of T = ^fim p v% ir / 'ks ~ 2 x 
10 4 /io.6(M/10 8 M Q ) 2 / 3 [(l + z)/10] K, where /jt = 0.6^0.6 
is a mean molecular mass of post-shock gas, and M is 



There is a typo in e-print version of paper 69] arXiv:1107.4094v2 
in formula (1): the right-hand side expression for T 2 (k) must be 
squared (private communication with authors) 
For the WDM initial power spectrum (1) from larXiv: 1 107.4094 v2 
(the square degree in the right-hand side of expression 1871 is 
omitted here) the quantitative accordance is essentially better 
(dashed lines). 
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the mass of halo or subhalo progenitors. The temper- 
ature of baryon matter gradually decreases afterwards 
due to the cooling processes (see @)- It results in frag- 
mentation to smaller clumps with Bonnor-Ebert mass, 
M BE ^ 700A/ Q (T/200K) 3 / 2 (n 6 /10 4 cm- 3 )- 1 /2 j where n& 

is the total number density of baryon particles. At the 
final stage of this fragmentation the stars and galaxies 
are formed (see Q and [J] for details). Since formation 
of galaxies is driven by gravity of dark matter, the spa- 
tial distribution of galaxies tracks the spatial distribu- 
tion of dark matter. In other words, the fluctuations of 
dark matter density 6dm (r) — Pdm(?)I~Pdm ~ 1 corre- 
late with fluctuations of galaxy number density 6 g {f) = 
n g {r)/n g — 1. The galaxy number density fluctuations 
have following Fourier amplitude in halo model, 



is appropriate 



8 iE 



1 



(N g \m) n(m)6z(m)y q (k, m)dm (88) 



(see Appendix IBl for details ). Here (N g \m) is a mean 
number of galaxies in halo with mass m, y g (k,m) is a 
Fourier-transform of galaxy number density profile and 
m m i n denotes the lowest limit for halo mass below which 
no galaxies are formed. Such limit naturally comes from 
degrading efficiency of star formation in halos of low mass 
and conditions imposed on the sample of galaxies (see 

HI). 

The considerations of previous subsection are summa- 
rized in following the galaxy power spectrum: 



p»w = 



1 



(Ng 



(27T) 3 j V , tg 

x [\y g (m,k)\ 2 ~W 2 (kR g )]dm 



n(m) 



(89) 



(N g \m) , w ._ . . 
— £ bi(m)n(m)y (m, k)dm 



Vl in {k), 



where R g = (3 (N g \m) /(47rn s )) 1 / 3 . The resulting equa- 
tion is similar to the corresponding expression for the 
galaxy power spectrum from [22j. The difference is 
caused by elimination of the noise as described above. 
Also the term (N„\m) has been obtained instead of 
(N g (N g — l)|m) in 22]. For large-mass halos the term 
(N g \m) is large, and it seems appropriate to assume the 
probability distribution p(N g \m) to be one of Poisson. 
In this case (N g \m) « (N g (N g — l)|m). However, this 
equation is not valid for low-mass halos. 

As it has been outlined in |Ty], galaxies within halo 
usually are disposed around the center (central galaxy) 
and within each of its subhalos (satellite galaxies). This 
gives a clue to find out how galaxies are distributed over 
the halo and how it is connected to the substructure. 
Massive halos usually undergo the violent relaxation. In 
such case, as follows from Eq. (fl3"l) . for the number den- 
sity of the satellite galaxies in halo the following equation 



i( r ) = ^2n s h(m,r) = y^ n (m)exp < 
3$(r) ' 



m$>{r) 



in>m' 



n„ exp 



,0 



kT 



-Mr), 

Ps 



(90) 



where n s h{m,r) denotes the dependence of the number 
density of halo particles (subhalos) with mass m on the 
radial distance. 

As above, we assume that galaxies are formed within 
subhalos with masses m > m' min , where m' min is less than 
m m i n because subhalos usually lose the mass due to tidal 
deprivation of their outskirts. Baryon matter (stars) are 
concentrated to the center and more tightly bound than 
dark matter and therefore is slightly changing, meanwhile 
dark matter is stripped off. Thus, the subhalo mass at 
the time of observation is probably a poor tracer for the 
potential well which determines galaxy properties such 
as stellar mass or luminosity. A better tracer is the sub- 
halo mass at the time when it falls into the host halo, 
or its maximum mass over its history. For massive ha- 
los with numerous satellites the presence of the central 
galaxy can be neglected. In such case, as follows from 
(f9"TJ)) . the assumption y g (m, k) ~ y(m, k) is correct. This 
result agrees with [73j, the spatial distribution of satel- 
lite galaxies using SDSS spectroscopic and photometric 
galaxy catalogs is studied therein. They found that satel- 
lite profiles generally have a universal form that well- 
fitted by NFW approximation. 

However, as long as low-mass halos possess a small 
number of galaxies, slow relaxation is also important for 
them. In result, the profile of satellite galaxy number 
density generally deviates from the profile of dark matter 
density. However, such discrepancy is difficult to detect 
because of large statistical errors in determination of the 
profile of the galaxy number density in such halos. Let us 
note that in this case we also can not discard the presence 
of central galaxy. 

The spatial number density of the galaxies is a sum of 
the halo and subhalo number densities: n g (f) — nhif) + 
n s h (r) . The spatial fluctuation of galaxy number density 
can be thereby split into fluctuations of halo and subhalo 
numbers densities: 

= n h (r) + n sh (r)_ i 

Uh + Ush 

r nh6h(r) +n s h6 s h(r)], (91) 



n h + n sh 



where as before the overlines denote averaging in space. 
Corresponding Fourier amplitude takes the form 



*-,£ = 



1 



9 l fc n h +n s h 



n hk\k+ n sh6 sh \j: 



(92) 



nh + n sh 



[1 + (N sh \m) y sh {k, m)] n(rn)6Am)dm, 
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FIG. 12: The halo model galaxy (solid line), dark matter (dotted line) power spectra calculated for ACDM cosmology with 
(f2 m , f^A, h, as, Us) = (0.26, 0.74, 0.72, 0.77, 0.95). The squares and triangles represent the observed galaxy power spectrum 
drawn from PSCz 25] and SDSS 26] galaxy catalogs respectively. 



where (N s h\m) is the average number of subhalos con- 
fined within the halo of mass m, it equals the num- 
ber of satellites. Since the average number of galax- 
ies accounts for central galaxy and satellites, therefore 
(N g \m) = 1 + (N sh \m). By comparing Eqs. ([55) and 

HMD, 



= nh + n s h= /[I + (N sh \m)]n(m)dm, (93) 



and 



y q (k,m) 



(N sh \m) y sh (k, m) + yjk, m) 



(N 8h \m) 



1 



(94) 



where y sh {k, m) is a Fourier image of the subhalo nunrber 
density profile, meanwhile y c (k,m) is a Fourier image of 
the probability of finding the central galaxy within the 
halo. 

If the central galaxies in all halos with masses m are 
located strictly in their centers, then y c (k,m) = 1. For 
massive halo, (N s h\m) ^> 1, so y g (k,m) ~ y sh (k,m), 
whereas for low-mass halo, (N s h\m) <C 1, y Q (k,m) ~ 
y c (k,m). As follows from Eq. (|90| for massive halos we 
can assume y sh (k,m) ~ y(k,m). For simplicity, let us 
extend this approach to the case of low-mass halo. It 
should not bring significant errors to y(k,m) because in 
the case of (N s h\m) <C 1 the central galaxy is dominating, 
so y g (k,m) ~y c (k,m). 



To specify the dependence (N s h\m) we can use the sub- 
halo mass function (|75p . In such case the number of 
satellite galaxies in halo with mass m is 



(N sh \m)= l-r-rdm' 
dm 



(95) 



the masses m and m' min are to be determined within the 
virial density regions with A vm p cr . 

To provide a direct link to the galaxy sample the 
CLF (IM! or CMF 3) can be also used instead 
of the subhalo mass function. The CLF, <b(L\m)dL, 
yields the average number of galaxies with luminosity 
L ± dL/2 which reside within a halo of mass m. The 
CMF, <5>{m*\m)dm*, yields the average number of galax- 
ies with stellar masses in the range m* ± dm*/2 which 
reside within halo of mass m. The CMF (as well as 
CLF) can be split into central and satellite parts so that 
$(m*|m) = ^(m, \m) + <r> c (m»|m). This allows us to 
calculate the average number of satellites with a stellar 
masses exceeding m* within the halo with mass m (see 
|16{ for details), 



(N sh \m,m*) = / $ s (m'*|m)dm'*, 



and the probability of finding the appropriate central 
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galaxy is 



(N c \m,m* 



$c(rn'Jm)dm^ 



where upper limit is assigned with infinity, although it 
actually does not exceed the halo mass to. To calcu- 
late the power spectrum of galaxies, we assume that 
(N g \m) = (N sh \m,m*) + (N c \m,m*) and 



_ (N sh \m,m*)y sh (k,m) + (N c \m,m*) y c (k,m) 

VgK ,m) (N sh \m,m*) + (N c \m,m*) 

(96) 
The average number of galaxies with a stellar mass 
larger than m* is given by 



[(N s h\m,m*) + (N c \m,m*}]n(m)dm. (97) 



The same calculations are valid for CLF. Hence, the halo 
model describes the connection between galaxy power 
spectrum and stellar masses or luminosities of the sample 
of galaxies. 

For example, we have calculated the galaxy power 
spectrum by using of CMF from [Ig. The CMF was 
used with scatter because it provides better consistency 
with observations. It was evaluated for ACDM cosmology 
with parameters (fi m , Qa, h, as, n s ) — (0.26, 0.74, 0.72, 
0.77, 0.95), so the same parameters were used herein. 
The CAMB code [zl, [75| was used to calculate the ini- 
tial dark matter power spectrum Vu n {k), with fif, = 0.05. 
For the galaxy power spectrum evaluation the Eq. (|9"6")l 
and Eq. (|59"j) were used with W(kR g ) replaced by 

lim v„(k, to) = 

{N sh \m,m*)W(kR s ) + (N c \m,m*)W(kR c ) 
(N sh \m,m*) + (N c \m,m*) 

where R s = (3 (N sh \m,m*) /(47m g )) 1/3 and R c = 
(3(W c |m,m*>/(4^n 9 )) 1 / 3 . 
We assumed 



y s h( k > m ) 



y(k,r s ,c')p(c'\m, z)dc' . 



allc 



where y(k,r s ,c) denotes the dependence (|84|) and 
p(c\m, z) is the function of concentration distribution 
(|4"9"]l with variance a c = 0.25. 

The calculated galaxy and dark matter power spec- 
tra are presented in Fig. 112] along with observed galaxy 
power spectra obtained from PSCz [25] and SDSS [26] 
galaxy catalogs. 

The upper solid line represents the assumption that 
the central galaxies in all halos with masses to are lo- 
cated strictly in their centers, so y c {k, m) — 1. The lower 
solid line represents the result for assumption that central 



galaxies are homogeneously distributed over the spheri- 
cal volume of radius l.lr s , so y c (k,m) = W(l.lr s k). We 
define the lower limit on the stellar masses of the galaxies 
to be m« = 5 x 10 6 M©. 

Thus, at large scales, k < 1 h _1 Mpc, the dark mat- 
ter and galaxy power spectra coincide, at galaxy cluster 
scales, 1 < k < 20 h _1 Mpc, they are close and start to di- 
verge at smaller scales, k > 20 h _1 Mpc, where luminous 
matter is substantially more clustered than dark matter. 



V. CONCLUSIONS 

The presented semi-analytical treatment is our imple- 
mentation of halo model capable to describe and interpret 
the clustering of the matter at the non-linear stage of evo- 
lution, both in simulated and observed Universe. Some 
of basic elements of theory is reviewed and improved to 
calculate the dark matter and galaxy power spectrum. 

At first we have analyzed the evolution of spherical 
scalar perturbations in different cosmologies for arbitrary 
initial density profile to find the relations between am- 
plitude of linear perturbations and overdensities at non- 
linear stages. These relations enable to estimate the mo- 
ments of turnaround, collapse and virialization of spher- 
ical overdensity for any set of main cosmological param- 
eters, f2 m , 17a, &k and Hq. The minimal amplitudes are 
estimated of initial density perturbations 5 m i n which can 
seed virialized halo in fl^ > or/and fl\ > cosmolo- 
gies. The exact analytical expressions (|3U|) . (|3"2")l , (|3"5j) . 
(J4H)) are suggested as well as analytical approximations 
([33]), dSZ]), g2). 

We propose a new technique for computation of halo 
concentration parameter, c, with phenomenology of halo 
merging, density profiles and statistical properties taken 
into account. The simple expression for estimation (|66[) 
depends on the relation of the halo overdensity, A vc or 
A„ c , and corresponding characteristic halo overdensity, 
A sc or A sc respectively. This relation is evaluated with- 
out computing redshift of halo collapse, z co i, by set of 
equations: (|5B |l . ([61 ]) . (|63"]l . (|30 ]l and ([BT jl or (p jl as well. 
Such technique has been applied to calculate the concen- 
tration parameter for ACDM and AWDM cosmological 
models, and the concordance with data of simulations 
[44L 1-451 ] for vast range of halo masses (Figs. El [7J) has 
been revealed. 

The parameters of Sheth-Tormen approximation for 
halo mass function were re-evaluated as p = 0.32 and 
q = 0.76 (see Fig. [8]) to provide best-fit to the data 
of GIF/ Virgo N-body simulations [59} (see Fig. [5]). 
The proposed approximations for subhalo mass function, 
Eq. ([75]) match the derived from GIF2 [27| (Fig. B). It 
implies that halo mass consists only by 30% of subha- 
los and the rest 70% is a diffuse dark matter which not 
to be associated with any substructure. However, we 
stress that resolving of substructure is limited by a some 
specified resolution level caused by finite resolution of 
numerical simulations. 
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Also, this modified and extended halo model enables 
to predict the dark matter and galaxy power spectra at 
small scales up to k ~ 100 h _1 Mpc by means of semi- 
analytical methods Eqs. ([8"5]). (j8"6"1) . (j8"9l . The estimated 
spectra agree with with non-linear power spectra deter- 
mined from Large Box and GIF2 N-body simulations 
(Fig. ITOl) as well as with estimations by galaxy catalogs 
PSCz [H and SDSS [26| (Fig. HU). Moreover, with the 
assumption on presence of the central galaxies in all ha- 
los with masses to (y c (k,m) — 1) the technique predicts 
galaxy power spectrum, matching well the observational 
one up to k ~ 20 h _1 Mpc. Meanwhile, when non-central 
position of most massive galaxies in halos is assumed, 
(y c (k,m) — W(l.lr a k)), the predictions agree with the 
observational data up to k ~ 80 h _1 Mpc. 

The calculated non-linear galaxy power spectrum for 
ACDM cosmology with (fi TO , fi A , h, a 8 , n s ) = (0.26, 0.74, 
0.72, 0.77, 0.95) corresponds to the observational one 
for lower limitation on the stellar masses of the galaxies 
to* = 5 • 10 6 Mq. To attain the same level of agreement 
of the predicted galaxy power spectrum with extracted 
from galaxy surveys at smaller scales (k > 80 h _1 Mpc), 
a new complicated approach for the formation of groups 
of galaxies must be elaborated. 

Despite the ambiguities in the definition of halo, de- 
termining of their mass, concentration and substructure, 
halo model provides a good reproduction of such char- 
acteristics of large-scale structure of the Universe as the 



power spectrum and correlation function of the spatial 
distribution of dark matter and galaxies. In this paper 
we show how communication between statistics of the 
dark matter clustering obtained from numerical simula- 
tions and galaxy statistics obtained from large galaxy 
surveys allows to calculate the power spectrum of the 
spatial distribution of the galaxies. 
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Appendix A: Fourier modes of dark matter density 
inhomogeneities 

The power spectrum can be derived in more rigorous 
manner by the series of following mathematical transfor- 
mations: 
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Here, the integration over the whole volume V has been 
split into integrations over volumes Vi, each occupied by 
spatially separated halos; the Fourier-transform of i-th 
density profile we denote as yt(k), it is normalized by 
their masses vrn. The halos are binned into the subsets 
with equal masses rrij , the number of halos is denoted by 
Nj. Also, it was assumed, that the halos of equal masses 
have identical density profiles, and correspondingly, their 
Fourier-transforms, y(m, k), are identical too. The Sr(m) 



is notation of Fourier-amplitude of spatial distribution of 
halo with masses m. The summation has been changed 
to the integration. 

Appendix B: Fourier modes of galaxy number 
density inhomogeneities 



The Fourier amplitude for relative fluctuations of 
galaxy concentration takes the following form: 
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Here, as in Appendix the integration over whole 
volume V was divided into integrations over the number 
of volumes VJ, filled by spatially separated halos. The 
Fourier-images of profiles of concentration of galaxies in 
halo are normalized by their number N g u] the Fourier- 
image of i-th profile of galaxy concentration is denoted 
by y g \i{k). The halo was partitioned into subsets of equal 
masses m and normalized by means of index j m . 

It was assumed that halo of equal masses have iden- 
tical profiles of galaxy concentration, so their Fourier- 
images y Jm,k) are identical; The sets of halos with 
equal masses m were partitioned into subsets containing 
the same number N of galaxies and denoted as n(m, N). 
The halo concentration n(m, N) is represented as a prod- 



uct of concentration of all halos with mass m and con- 
ditional probability of event that these halos contain N 
galaxies each, n(m,N) — n(m)p(N\m). The designation 
was introduced for Fourier-amplitude of spatial distribu- 
tion of halos with masses m and containing N galax- 
ies as S^(m,N). It was assumed, that spatial distri- 
bution of halos of mass m and number of galaxies N 
match the spatial distribution of all halos with masses 
m, St(m) = Sj;(m,N). The average number of galaxies 

oo 

in halo of mass m is denoted as (N\m) = Yl Np(N\m), 

N=0 

also the sum was replaced by intergration. 
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